library(leaflet)
library(htmlwidgets)
library(plotly)
## Caricamento del pacchetto richiesto: ggplot2
##
## Caricamento pacchetto: 'plotly'
## Il seguente oggetto è mascherato da 'package:ggplot2':
##
## last_plot
## Il seguente oggetto è mascherato da 'package:stats':
##
## filter
## Il seguente oggetto è mascherato da 'package:graphics':
##
## layout
library(leaflet.extras)
library(dplyr)
##
## Caricamento pacchetto: 'dplyr'
## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag
## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, union
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(htmltools) # Assicurati di caricare il pacchetto htmltools
library(leafpop)
library(gridExtra)
##
## Caricamento pacchetto: 'gridExtra'
## Il seguente oggetto è mascherato da 'package:dplyr':
##
## combine
library(stats)
# Carica il nuovo file combinato
NDVI_VALUES <- st_read("G:/Altri computer/Il_mio_computer/DOTTORATO/PROGETTI/OLIVASTRO_PAULILATINO/REGRESSIONE/VETTORIALI/NDVI_VALUES/merged_data_ndvi.shp")
## Reading layer `merged_data_ndvi' from data source
## `G:\Altri computer\Il_mio_computer\DOTTORATO\PROGETTI\OLIVASTRO_PAULILATINO\REGRESSIONE\VETTORIALI\NDVI_VALUES\merged_data_ndvi.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 60928 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 473672.7 ymin: 4433508 xmax: 481731 ymax: 4437127
## Projected CRS: WGS 84 / UTM zone 32N
# Imposta l'opzione scipen su un valore elevato per eliminare la notazione esponenziale dei valori
options(scipen = 999)
# Carico il nuovo file con le classi
pixel_trend <- st_read("G:/Altri computer/Il_mio_computer/DOTTORATO/PROGETTI/OLIVASTRO_PAULILATINO/REGRESSIONE/VETTORIALI/NDVI_VALUES/pixels_trend2.shp")
## Reading layer `pixels_trend2' from data source
## `G:\Altri computer\Il_mio_computer\DOTTORATO\PROGETTI\OLIVASTRO_PAULILATINO\REGRESSIONE\VETTORIALI\NDVI_VALUES\pixels_trend2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 859 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 8.691293 ymin: 40.05143 xmax: 8.785804 ymax: 40.08405
## Geodetic CRS: WGS 84
# Carico il file della geolocalizzazione dei campioni
samples <- st_read("G:/Altri computer/Il_mio_computer/DOTTORATO/PROGETTI/OLIVASTRO_PAULILATINO/REGRESSIONE/VETTORIALI/sample_points.shp")
## Reading layer `sample_points' from data source
## `G:\Altri computer\Il_mio_computer\DOTTORATO\PROGETTI\OLIVASTRO_PAULILATINO\REGRESSIONE\VETTORIALI\sample_points.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 86 features and 12 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: 8.690753 ymin: 40.02248 xmax: 8.8019 ymax: 40.11557
## Geodetic CRS: WGS 84
lim_paul_wgs84 <- st_read("G:/Altri computer/Il_mio_computer/DOTTORATO/PROGETTI/OLIVASTRO_PAULILATINO/VETTORIALI/Limite_Amministrativo_Paulilatino_wgs84.shp")
## Reading layer `Limite_Amministrativo_Paulilatino_wgs84' from data source
## `G:\Altri computer\Il_mio_computer\DOTTORATO\PROGETTI\OLIVASTRO_PAULILATINO\VETTORIALI\Limite_Amministrativo_Paulilatino_wgs84.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 8.683562 ymin: 40.0118 xmax: 8.815261 ymax: 40.1325
## Geodetic CRS: WGS 84
focolai_wgs84 <- st_read("G:/Altri computer/Il_mio_computer/DOTTORATO/PROGETTI/OLIVASTRO_PAULILATINO/REGRESSIONE/VETTORIALI/FOCOLAI_wgs84.shp")
## Reading layer `FOCOLAI_wgs84' from data source
## `G:\Altri computer\Il_mio_computer\DOTTORATO\PROGETTI\OLIVASTRO_PAULILATINO\REGRESSIONE\VETTORIALI\FOCOLAI_wgs84.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 8.691165 ymin: 40.00848 xmax: 8.801279 ymax: 40.07877
## Geodetic CRS: WGS 84
# Define a custom color palette for classes
class_palette <- c("Crescente" = "#00ff00", # Green
"Stazionario" = "#adff2f", # Yellow-green
"Decrescente Lieve" = "yellow",
"Decrescente Moderato" = "orange",
"Decrescente Marcato" = "red")
# Create a color factor with the custom palette and class labels
color_factor <- leaflet::colorFactor(palette = class_palette, levels = c("Crescente", "Stazionario", "Decrescente Lieve", "Decrescente Moderato", "Decrescente Marcato"))
# Definisci l'ordine desiderato delle etichette delle classi
custom_order <- c("Crescente", "Stazionario", "Decrescente Lieve", "Decrescente Moderato", "Decrescente Marcato")
# Assicurati che 'classe_trend' sia un factor con l'ordine definito
pixel_trend$clss_tr <- factor(pixel_trend$clss_tr, levels = custom_order)
# Crea una funzione di colorazione per i cerchi
color_factor_circle <- colorFactor(
palette = c("green", "red"),
domain = c("+", "-")
)
# Ottieni i valori unici di COD
cod_values <- unique(NDVI_VALUES$COD)
# Funzione per creare il codice HTML del grafico
crea_grafico_html <- function(cod) {
# Filtra i dati per il valore corrente di COD
dati <- NDVI_VALUES[NDVI_VALUES$COD == cod, ]
# Crea la serie temporale
serie_temporale <- ts(dati$ndvi, start = c(2018, 1), end = c(2023, 8), frequency = 12)
# Crea il grafico con plot.ts
png(filename="temp.png")
plot.ts(serie_temporale, main=cod, ylab='NDVI')
lines(smooth.spline(time(serie_temporale), serie_temporale)$y, col = "blue") # n/100 for seasonal component
lines(supsmu(time(serie_temporale), serie_temporale, span=0.5)$y, col = "red") # n/2 for trend component
dev.off()
# Leggi l'immagine e codificala in base64
img <- base64enc::dataURI(file = "temp.png", mime = "image/png")
# Rimuovi il file temporaneo
file.remove("temp.png")
# Ritorna il codice HTML dell'immagine
return(paste0('<img src="', img, '"/>'))
}
# Crea la mappa Leaflet
map <- leaflet(options = leafletOptions(maxZoom = 22)) %>%
addProviderTiles("Esri.WorldImagery", group = "Imagery", options = providerTileOptions()) %>%
addPolygons(data = lim_paul_wgs84,
fillOpacity = 0,
color = "black",
weight = 2,
group = 'Lim. Amm. Paulilatino' # Add a group name
) %>%
addPolygons(data = focolai_wgs84,
fillOpacity = 0,
color = "red",
weight = 2,
group = 'Outbreaks'
) %>%
addCircleMarkers(data = samples,
lng = ~long,
lat = ~lat,
fillColor = ~color_factor_circle(samples$positivo),
fillOpacity = 0.6,
popup = ~paste("AREA:", location, "<br/>CAMPIONE:", id_sample, "<br/>SINTOMATICO:", sin_asin, "<br/>POSITIVITÀ:", positivo),
radius = 3,
group = 'samples',
stroke = FALSE
) %>%
addLegend(pal = color_factor_circle,
values = samples$positivo,
title = "Sample Positivity",
opacity = 0.6,
position = "bottomright"
) %>%
addPolygons(data = pixel_trend,
fillColor = ~color_factor(clss_tr),
color = "black",
fillOpacity = 0.6,
highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE),
popup =~paste(
"COD:", COD, "<br/>coef:", coef, "<br/>classe:", clss_tr),
group ='Crowns',
stroke=TRUE,
weight=1
) %>%
addLegend(title="Trend: lm(NDVI ~ Month)",
pal=color_factor,
values=pixel_trend$clss_tr,
opacity=0.6,
position="topright"
) %>%
addFullscreenControl() %>%
addLayersControl(
overlayGroups=c("Crowns","samples","Lim. Amm. Paulilatino","Outbreaks"),
options=layersControlOptions(collapsed=TRUE)
)
# Add the geolocation control
map <- activateGPS(map) %>%
addControlGPS(
options=gpsOptions(
position="topleft",
autoCenter=TRUE,
))
# Aggiungi i popup con i grafici alla mappa Leaflet
for (i in seq_along(cod_values)) {
cod <- cod_values[i]
# Filtra i dati di NDVI_VALUES per il valore corrente di COD
NDVI_VALUES_cod <- NDVI_VALUES[NDVI_VALUES$COD == cod, ]
lng <- st_coordinates(NDVI_VALUES_cod$geometry)[, 1]
lat <- st_coordinates(NDVI_VALUES_cod$geometry)[, 2]
# Crea il popup con il grafico HTML
html <- crea_grafico_html(cod)
# Aggiungi il popup al marker sulla mappa
map <- addMarkers(map, lng = lng, lat = lat, popup = html)
}
map